home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / cagd_lib / cagdalen.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  10.5 KB  |  357 lines
  1. /******************************************************************************
  2. * CagdALen.c - arc length functions and unit length normalization.          *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Apr. 93.                          *
  5. ******************************************************************************/
  6.  
  7. #include "cagd_loc.h"
  8.  
  9. #define MAX_ALEN_IMPROVE_ITERS 5
  10. #define ARCLEN_CONST_SET_EPSILON 0.001
  11.  
  12. /******************************************************************************
  13. * Subdivides the given curves to curves, each with less of control polygon    *
  14. * less than or equal to MaxLen.    Returned is a list of curves.              *
  15. ******************************************************************************/
  16. CagdCrvStruct *CagdLimitCrvArcLen(CagdCrvStruct *Crv, CagdRType MaxLen)
  17. {
  18.     if (CagdCrvArcLenPoly(Crv) > MaxLen) {
  19.     CagdRType TMin, TMax;
  20.     CagdCrvStruct *Crv1, *Crv2, *Crv1MaxLen, *Crv2MaxLen;
  21.  
  22.     CagdCrvDomain(Crv, &TMin, &TMax);
  23.  
  24.     Crv1 = CagdCrvSubdivAtParam(Crv, (TMin + TMax) / 2.0);
  25.     Crv2 = Crv1 -> Pnext;
  26.  
  27.     Crv1MaxLen = CagdLimitCrvArcLen(Crv1, MaxLen);
  28.     Crv2MaxLen = CagdLimitCrvArcLen(Crv2, MaxLen);
  29.  
  30.     CagdCrvFree(Crv1);
  31.     CagdCrvFree(Crv2);
  32.  
  33.     for (Crv1 = Crv1MaxLen; Crv1 -> Pnext != NULL; Crv1 = Crv1 -> Pnext);
  34.     Crv1 -> Pnext = Crv2MaxLen;
  35.     return Crv1MaxLen;
  36.     }
  37.     else
  38.         return CagdCrvCopy(Crv);
  39. }
  40.  
  41. /******************************************************************************
  42. * Computes a bound on the arc length of a curve by computing the length of    *
  43. * its control polygon.                                  *
  44. ******************************************************************************/
  45. CagdRType CagdCrvArcLenPoly(CagdCrvStruct *Crv)
  46. {
  47.     int i;
  48.     CagdCrvStruct
  49.     *CrvE3 = CagdCoerceCrvTo(Crv, CAGD_PT_E3_TYPE);
  50.     CagdRType
  51.     Len = 0.0,
  52.     **Points = CrvE3 -> Points;    
  53.  
  54.     for (i = 1; i < CrvE3 -> Length; i++)
  55.         Len += sqrt(SQR(Points[1][i] - Points[1][i - 1]) +
  56.             SQR(Points[2][i] - Points[2][i - 1]) +
  57.             SQR(Points[3][i] - Points[3][i - 1]));
  58.  
  59.     CagdCrvFree(CrvE3);
  60.  
  61.     return Len;
  62. }
  63.  
  64. /******************************************************************************
  65. * Normalizes the given curve to be a unit length curve, by computing a scalar *
  66. * curve to multiply with this curve. Returns the multiplied curve if Mult, or *
  67. * just the scalar curve, otherwise.                          *
  68. ******************************************************************************/
  69. CagdCrvStruct *CagdCrvUnitLenScalar(CagdCrvStruct *OrigCrv, CagdBType Mult,
  70.                                CagdRType Epsilon)
  71. {
  72.     int i, j;
  73.     CagdCrvStruct
  74.     *ScalarCrv = NULL,
  75.     *Crv = CAGD_IS_BEZIER_CRV(OrigCrv) ? CnvrtBezier2BsplineCrv(OrigCrv)
  76.                        : CagdCrvCopy(OrigCrv);
  77.     CagdBType
  78.         IsRational = CAGD_IS_RATIONAL_CRV(Crv);
  79.  
  80.     for (i = 0; i < MAX_ALEN_IMPROVE_ITERS; i++) {
  81.     CagdCrvStruct *UnitCrvAprx, *ScalarCrvSqr,
  82.         *DotProdCrv = CagdCrvDotProd(Crv, Crv);
  83.     CagdRType Min, Max, *SCPoints,
  84.         *DPPoints = DotProdCrv -> Points[1];
  85.  
  86.     if (ScalarCrv)
  87.         CagdCrvFree(ScalarCrv);
  88.     ScalarCrv = CagdCrvCopy(DotProdCrv);
  89.     SCPoints = ScalarCrv -> Points[1];
  90.  
  91.     /* Compute an approximation to the inverse function. */
  92.     for (j = 0; j < ScalarCrv -> Length; j++, DPPoints++) {
  93.         *SCPoints++ = *DPPoints > 0.0 ? 1.0 / sqrt(*DPPoints) : 1.0;
  94.     }
  95.     ScalarCrvSqr = CagdCrvMult(ScalarCrv, ScalarCrv);
  96.  
  97.     /* Multiply the two to test the error. */
  98.     UnitCrvAprx = CagdCrvMult(ScalarCrvSqr, DotProdCrv);
  99.     CagdCrvFree(ScalarCrvSqr);
  100.  
  101.     CagdCrvMinMax(UnitCrvAprx, 1, &Min, &Max);
  102.  
  103.     if (1.0 - Min < Epsilon && Max - 1.0 < Epsilon) {
  104.         CagdCrvFree(UnitCrvAprx);
  105.         CagdCrvFree(DotProdCrv);
  106.         break;
  107.     }
  108.     else {
  109.         /* Refine in regions where the error is too large. */
  110.         int k,
  111.             Length = UnitCrvAprx -> Length,
  112.             Order = UnitCrvAprx -> Order,
  113.             KVLen = Length + Order;
  114.         CagdRType
  115.         *KV = UnitCrvAprx -> KnotVector,
  116.         *RefKV = IritMalloc(sizeof(CagdRType) * 2 * Length),
  117.             *Nodes = BspKnotNodes(KV, KVLen, Order),
  118.         **Points = UnitCrvAprx -> Points;
  119.  
  120.         for (j = k = 0; j < Length; j++) {
  121.         CagdRType
  122.             V = 1.0 - (IsRational ? Points[1][j] / Points[0][j]
  123.                       : Points[1][j]);
  124.  
  125.         if (ABS(V) > Epsilon) {
  126.             int Index = BspKnotLastIndexLE(KV, KVLen, Nodes[j]);
  127.  
  128.             if (APX_EQ(KV[Index], Nodes[j])) {
  129.             if (j > 0)
  130.                 RefKV[k++] = (Nodes[j] + Nodes[j - 1]) / 2.0;
  131.             if (j < Length - 1)
  132.                 RefKV[k++] = (Nodes[j] + Nodes[j + 1]) / 2.0;
  133.             }
  134.             else
  135.             RefKV[k++] = Nodes[j];
  136.         }
  137.         }
  138.         CagdCrvFree(UnitCrvAprx);
  139.         CagdCrvFree(DotProdCrv);
  140.  
  141.         IritFree((VoidPtr) Nodes);
  142.  
  143.         if (k == 0) {
  144.         /* No refinement was found needed - return current curve. */
  145.         IritFree((VoidPtr) RefKV);
  146.         break;
  147.         }
  148.         else {
  149.         CagdCrvStruct
  150.             *CTmp = CagdCrvRefineAtParams(Crv, FALSE, RefKV, k);
  151.  
  152.         IritFree((VoidPtr) RefKV);
  153.         CagdCrvFree(Crv);
  154.         Crv = CTmp;
  155.         }
  156.  
  157.     }
  158.     }
  159.  
  160.     CagdCrvFree(Crv);
  161.  
  162.     /* Error is probably within bounds - returns this unit length approx. */
  163.     if (Mult) {
  164.     CagdCrvStruct *CrvX, *CrvY, *CrvZ, *CrvW, *CTmp;
  165.     int MaxCoord = CAGD_NUM_OF_PT_COORD(OrigCrv -> PType);
  166.  
  167.     CagdCrvSplitScalar(ScalarCrv, &CrvW, &CrvX, &CrvY, &CrvZ);
  168.     CagdCrvFree(ScalarCrv);
  169.     ScalarCrv = CagdCrvMergeScalar(CrvW,
  170.                        CrvX,
  171.                        MaxCoord > 1 ? CrvX : NULL,
  172.                        MaxCoord > 2 ? CrvX : NULL);
  173.     CagdCrvFree(CrvX);
  174.     if (CrvW)
  175.         CagdCrvFree(CrvW);
  176.  
  177.     CTmp = CagdCrvMult(ScalarCrv, OrigCrv);
  178.     CagdCrvFree(ScalarCrv);
  179.  
  180.     return CTmp;
  181.     }
  182.     else {
  183.     return ScalarCrv;
  184.     }
  185. }
  186.  
  187. /******************************************************************************
  188. * Computes the curve which is a square root approximation to a given scalar   *
  189. * curve, to within epsilon.                              *
  190. ******************************************************************************/
  191. CagdCrvStruct *CagdCrvSqrtScalar(CagdCrvStruct *OrigCrv, CagdRType Epsilon)
  192. {
  193.     int i, j;
  194.     CagdCrvStruct
  195.     *ScalarCrv = NULL,
  196.     *Crv = CAGD_IS_BEZIER_CRV(OrigCrv) ? CnvrtBezier2BsplineCrv(OrigCrv)
  197.                        : CagdCrvCopy(OrigCrv);
  198.     CagdBType
  199.         IsRational = CAGD_IS_RATIONAL_CRV(Crv);
  200.  
  201.     for (i = 0; i < MAX_ALEN_IMPROVE_ITERS; i++) {
  202.     CagdCrvStruct *ScalarCrvSqr, *ErrorCrv;
  203.     CagdRType Min, Max, *SCPoints,
  204.         *Points = Crv -> Points[1],
  205.         *WPoints = IsRational ? Crv -> Points[0] : NULL;
  206.  
  207.     if (ScalarCrv)
  208.         CagdCrvFree(ScalarCrv);
  209.     ScalarCrv = CagdCrvCopy(Crv);
  210.     SCPoints = ScalarCrv -> Points[1];
  211.  
  212.     /* Compute an approximation to the inverse function. */
  213.     for (j = 0; j < ScalarCrv -> Length; j++) {
  214.         CagdRType
  215.         V = IsRational ? *Points++ / *WPoints++ : *Points++;
  216.  
  217.         *SCPoints++ = V >= 0.0 ? sqrt(V) : 0.0;
  218.     }
  219.     ScalarCrvSqr = CagdCrvMult(ScalarCrv, ScalarCrv);
  220.  
  221.     /* Compare the two to test the error. */
  222.     ErrorCrv = CagdCrvSub(ScalarCrvSqr, Crv);
  223.     CagdCrvFree(ScalarCrvSqr);
  224.  
  225.     CagdCrvMinMax(ErrorCrv, 1, &Min, &Max);
  226.  
  227.     if (Min > -Epsilon && Max < Epsilon) {
  228.         CagdCrvFree(ErrorCrv);
  229.         break;
  230.     }
  231.     else {
  232.         /* Refine in regions where the error is too large. */
  233.         int k,
  234.             Length = ErrorCrv -> Length,
  235.             Order = ErrorCrv -> Order,
  236.             KVLen = Length + Order;
  237.         CagdRType
  238.         *KV = ErrorCrv -> KnotVector,
  239.         *RefKV = IritMalloc(sizeof(CagdRType) * 2 * Length),
  240.             *Nodes = BspKnotNodes(KV, KVLen, Order),
  241.         **ErrPoints = ErrorCrv -> Points;
  242.  
  243.         for (j = k = 0; j < Length; j++) {
  244.         CagdRType
  245.             V = 1.0 - (IsRational ? ErrPoints[1][j] / ErrPoints[0][j]
  246.                       : ErrPoints[1][j]);
  247.  
  248.         if (ABS(V) > Epsilon) {
  249.             int Index = BspKnotLastIndexLE(KV, KVLen, Nodes[j]);
  250.  
  251.             if (APX_EQ(KV[Index], Nodes[j])) {
  252.             if (j > 0)
  253.                 RefKV[k++] = (Nodes[j] + Nodes[j - 1]) / 2.0;
  254.             if (j < Length - 1)
  255.                 RefKV[k++] = (Nodes[j] + Nodes[j + 1]) / 2.0;
  256.             }
  257.             else
  258.             RefKV[k++] = Nodes[j];
  259.         }
  260.         }
  261.         CagdCrvFree(ErrorCrv);
  262.  
  263.         IritFree((VoidPtr) Nodes);
  264.  
  265.         if (k == 0) {
  266.         /* No refinement was found needed - return current curve. */
  267.         IritFree((VoidPtr) RefKV);
  268.         break;
  269.         }
  270.         else {
  271.         CagdCrvStruct
  272.             *CTmp = CagdCrvRefineAtParams(Crv, FALSE, RefKV, k);
  273.  
  274.         IritFree((VoidPtr) RefKV);
  275.         CagdCrvFree(Crv);
  276.         Crv = CTmp;
  277.         }
  278.  
  279.     }
  280.     }
  281.  
  282.     CagdCrvFree(Crv);
  283.  
  284.     return ScalarCrv;
  285. }
  286.  
  287. /******************************************************************************
  288. * Normalizes the given curve to be a unit length curve, by computing a scalar *
  289. * curve to multiply with this curve. Returns the composed curve if Compose,   *
  290. * or just the scalar curve, otherwise.                          *
  291. * Note that the returned curves is approximately constant speed curve and     *
  292. * not arc-length, since a Bezier curve always has a domain 0 to 1 nomatter    *
  293. * the size of the curve itself.                              *
  294. ******************************************************************************/
  295. CagdCrvStruct *CagdCrvArcLenCrv(CagdCrvStruct *Crv, CagdRType Epsilon)
  296. {
  297.     CagdCrvStruct
  298.     *DerivCrv = CagdCrvDerive(Crv),
  299.         *DerivMagSqrCrv = CagdCrvDotProd(DerivCrv, DerivCrv),
  300.     *DerivMagCrv = CagdCrvSqrtScalar(DerivMagSqrCrv, Epsilon),
  301.         *ArcLenCrv = CagdCrvIntegrate(DerivMagCrv);
  302.  
  303.     CagdCrvFree(DerivCrv);
  304.     CagdCrvFree(DerivMagSqrCrv);
  305.     CagdCrvFree(DerivMagCrv);
  306.  
  307.     return ArcLenCrv;
  308. }
  309.  
  310. /******************************************************************************
  311. * Computes a tight approximation to the arc length of a curve.              *
  312. ******************************************************************************/
  313. CagdRType CagdCrvArcLen(CagdCrvStruct *Crv, CagdRType Epsilon)
  314. {
  315.     CagdRType TMin, TMax, *Pt;
  316.     CagdCrvStruct
  317.     *ArcLenCrv = CagdCrvArcLenCrv(Crv, Epsilon);
  318.  
  319.     CagdCrvDomain(ArcLenCrv, &TMin, &TMax);
  320.     Pt = CagdCrvEval(ArcLenCrv, TMax);
  321.     CagdCrvFree(ArcLenCrv);
  322.  
  323.     return CAGD_IS_RATIONAL_CRV(ArcLenCrv) ? Pt[1] / Pt[0] : Pt[1];
  324. }
  325.  
  326. /******************************************************************************
  327. * Computes parameter values to move steps of Length at a time. Returned is a  *
  328. * list of parameter values to step at.                          *
  329. ******************************************************************************/
  330. CagdPtStruct *CagdCrvArcLenSteps(CagdCrvStruct *Crv, CagdRType Length,
  331.                                 CagdRType Epsilon)
  332. {
  333.     CagdRType TMin, TMax, *Pt, CrvLength, Len;
  334.     CagdPtStruct *Param,
  335.     *ParamList = NULL;
  336.     CagdCrvStruct
  337.     *ArcLenCrv = CagdCrvArcLenCrv(Crv, Epsilon);
  338.  
  339.     CagdCrvDomain(ArcLenCrv, &TMin, &TMax);
  340.     Pt = CagdCrvEval(ArcLenCrv, TMax);
  341.  
  342.     CrvLength = CAGD_IS_RATIONAL_CRV(ArcLenCrv) ? Pt[1] / Pt[0] : Pt[1];
  343.  
  344.     for (Len = CrvLength - Length; Len > 0.0; Len -= Length) {
  345.     Param = CagdCrvConstSet(ArcLenCrv, 1, ARCLEN_CONST_SET_EPSILON, Len);
  346.     if (Param == NULL || Param -> Pnext != NULL)
  347.         FATAL_ERROR(CAGD_ERR_REPARAM_NOT_MONOTONE);
  348.  
  349.     Param -> Pnext = ParamList;
  350.     ParamList = Param;
  351.     }
  352.  
  353.     CagdCrvFree(ArcLenCrv);
  354.  
  355.     return ParamList;
  356. }
  357.